/*********************************************************************
*                                                                    *
*  rkc -- integrates    dy/dt = f(y,t)     using runge-kutta         *
*                                                                    *
* Parameters:                                                        *
*      v_init_ptr     -- pointer to where the y initial is           *
*                        and where the y final goes                  *
*      h              -- time step                                   *
*      f              -- pointer to function which computes f(y,t)   *
*                                                                    *
* Returns: nothing                                                   *
*                                                                    *
* WARNING uses global variables for intermediate values              *
*********************************************************************/

#include "first_macro_2d.h"  /* Where global parameters are defined */
#include "structs_2d.h"      /* Where structures are defined */
#include "derivs_2d.h"       /* Where derivatives functions are defined */
#include "gen_2d.h"

void rkc(struct field_array *v_init_ptr,
       FLOAT h,  
       struct GRID_PAR *grid_ptr,
       struct FUNCTION_PAR *equation_par,
       void (* FF)(FF_ARG_DUMMY)) {

FLOAT x;      /* initial time */    
FLOAT xh;     /* initial time + hh */
FLOAT hh;     /* hh half the h */
FLOAT h6;     /* one sixth of h */
int n_fields = (*grid_ptr).n_fields;
long int n_var = ((*grid_ptr).final_grid_1-(*grid_ptr).start_grid_1)
               *((*grid_ptr).final_grid_2-(*grid_ptr).start_grid_2);

/****************** global variables defined above main *****************/


/* struct field_array dv_temp;  */      /* derivative temporal value         */ 
/* struct field_array dv_sum;   */      /* sum values of deriv               */ 
/* struct field_array v_temp;   */      /* intermediate value of v           */ 

/* struct field_array *dv_temp_ptr = &dv_temp;  */  /* derivative temporal value        */ 
/* struct field_array *dv_sum_ptr = &dv_sum;    */  /* sum          values of deriv     */ 
/* struct field_array *v_temp_ptr = &v_temp;    */  /* intermediate value of v          */ 


/******************* end global variables *******************************/


x = (*v_init_ptr).time;
hh = h*(0.5);
h6 = h/(6.0);
xh = x + hh;




 FF(grid_ptr,v_init_ptr,dv_temp_ptr,equation_par);       /* take derivatives at initial time */

/* put derivative in dev_sum = (k1 + 2k2 + 2k3 + k4) */
/* dev_sum = k1 */
#ifdef DEBUG
 printf("after first function evaluation in RKC \n");
fflush(stdout);
#endif

 {int field_ind;
 long int ind;
   for (field_ind = 0; field_ind < n_fields; field_ind++){ 
       for (ind = 0; ind < n_var; ++ind){
	   (*dv_sum_ptr).u[field_ind][ind] = (*dv_temp_ptr).u[field_ind][ind];
       }
   }
 }


/* asign values to v_temp (v_temp = v_init + hh*dv_init) */

 {int field_ind;
 long int ind;
   for (field_ind = 0; field_ind < n_fields; field_ind++){ 
       for (ind = 0; ind < n_var; ++ind){
   (*v_temp_ptr).u[field_ind][ind] = (*v_init_ptr).u[field_ind][ind] 
                          + (*dv_temp_ptr).u[field_ind][ind]*hh;
       }
   }
 }

/* asign time to v_temp */
(*v_temp_ptr).time = xh;

/* take again derivative at time xh = t + h/2 and put them in dv_temp */
 FF(grid_ptr,v_temp_ptr,dv_temp_ptr,equation_par);


/* put derivative in dev_sum = (k1 + 2k2 + 2k3 + k4) */
/* dev_sum = k1 + 2k2 */

 {int field_ind;
 long int ind;
   for (field_ind = 0; field_ind < n_fields; field_ind++){ 
       for (ind = 0; ind < n_var; ++ind){
	   (*dv_sum_ptr).u[field_ind][ind] = (*dv_sum_ptr).u[field_ind][ind] 
	       + 2*(*dv_temp_ptr).u[field_ind][ind];
       }
   }
 }



/* asign values to v_temp = v_init + (h/2)*dv_temp  */
 {int field_ind;
 long int ind;
   for (field_ind = 0; field_ind < n_fields; field_ind++){ 
       for (ind = 0; ind < n_var; ++ind){
	   (*v_temp_ptr).u[field_ind][ind] = (*v_init_ptr).u[field_ind][ind] 
	       + (*dv_temp_ptr).u[field_ind][ind]*hh;
       }
   }
 }

 /* time does not need to be updated */

 /* takes derivatives at x+h/2, v_temp (new) and put them in dv_med */
 FF(grid_ptr,v_temp_ptr,dv_temp_ptr,equation_par);


/* put derivative in dev_sum = (k1 + 2k2 + 2k3 + k4) */
/* dev_sum = k1 + 2k2 + 2k3 */

 {int field_ind;
 long int ind;
   for (field_ind = 0; field_ind < n_fields; field_ind++){ 
       for (ind = 0; ind < n_var; ++ind){
	   (*dv_sum_ptr).u[field_ind][ind] = (*dv_sum_ptr).u[field_ind][ind] 
	       +2*(*dv_temp_ptr).u[field_ind][ind];
       }
   }
 }


/* asign values to v_temp = v_init + h*dv_temp */
 {int field_ind;
 long int ind;   
   for (field_ind = 0; field_ind < n_fields; field_ind++){ 
     for (ind = 0; ind < n_var; ++ind){
	 (*v_temp_ptr).u[field_ind][ind] = (*v_init_ptr).u[field_ind][ind] 
	     + (*dv_temp_ptr).u[field_ind][ind]*h;
     }
   }
 }

 /* asigns new time to v_temp = x + h */
 (*v_temp_ptr).time = x + h;

/* takes derivatives at x+h with new v_temp and puts them in dv_temp */
 FF(grid_ptr,v_temp_ptr,dv_temp_ptr,equation_par);

/* asign value to v_fina = v_init + h6*(dv_temp + dv_temp)*/
 {int field_ind;
 long int ind;   
 for (field_ind = 0; field_ind < n_fields; field_ind++){ 
     for (ind = 0; ind < n_var; ++ind){
	 (*v_init_ptr).u[field_ind][ind] = (*v_init_ptr).u[field_ind][ind] 
	     + h6*((*dv_sum_ptr).u[field_ind][ind] 
		   +  (*dv_temp_ptr).u[field_ind][ind]);
     }
 }
 }

 /* asign returning time */
(*v_init_ptr).time = x + h;
/*  printf("time at RKC = %f",x+h); */
}
